library(Seurat)
library(princurve)
library(monocle)
library(gprofiler2)
library(seriation)
library(Matrix)
library(dplyr)
library(RColorBrewer)
library(ggplot2)
library(ggExtra)
library(cowplot)
library(wesanderson)
#Set ggplot theme as classic
theme_set(theme_classic())Hem.data <- readRDS("../QC.filtered.cells.RDS")
Idents(Hem.data) <- Hem.data$Cell_identDimPlot(object = Hem.data,
group.by = "Cell_ident",
reduction = "spring",
cols = c("#83c3b8", #"ChP"
"#009fda", #"ChP_progenitors"
"#68b041", #"Dorso-Medial_pallium"
"#e46b6b", #"Hem"
"#e3c148", #"Medial_pallium"
"#b7d174", #2
"grey40", #4
"black", #5
"#3e69ac" #"Thalamic_eminence"
))Neurons.data <- subset(Hem.data, idents = c("seurat_clusters_2"))
DimPlot(Neurons.data ,
reduction = "spring",
pt.size = 1,
cols = c("#b7d174")) + NoAxes()p1 <- FeaturePlot(object = Neurons.data ,
features = c("BP_signature1","LN_signature1"),
pt.size = 0.5,
cols = rev(brewer.pal(10,"Spectral")),
reduction = "spring",
order = T) & NoAxes()
p2 <- FeaturePlot(object = Neurons.data ,
features = c("Foxg1", "Trp73"),
pt.size = 0.5,
cols = c("grey90", brewer.pal(9,"YlGnBu")),
reduction = "spring",
order = T) & NoAxes()
p1 / p2Separation between the 2 lineage seems straightforward. We use louvain clustering to split the two.
Neurons.data <- RunPCA(Neurons.data, verbose = FALSE)
Neurons.data <- FindNeighbors(Neurons.data,
dims = 1:10,
k.param = 8)
Neurons.data <- FindClusters(Neurons.data, resolution = 0.05)## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
##
## Number of nodes: 2835
## Number of edges: 56608
##
## Running Louvain algorithm...
## Maximum modularity in 10 random starts: 0.9640
## Number of communities: 2
## Elapsed time: 0 seconds
DimPlot(Neurons.data,
reduction = "spring",
cols = c("#cc391b","#026c9a"),
pt.size = 0.5) & NoAxes()Neurons.data$Lineage <- sapply(as.numeric(Neurons.data$SCT_snn_res.0.05),
FUN = function(x) {x= c("Hem","Pallial")[x]})DimPlot(object = Neurons.data,
group.by = "Lineage",
reduction = "spring",
cols = c("#cc391b","#026c9a"),
pt.size = 0.5) & NoAxes()Trajectories.Hem <- Neurons.data@meta.data %>%
select("Barcodes", "nUMI", "Spring_1", "Spring_2", "Lineage") %>%
filter(Lineage == "Hem")fit <- principal_curve(as.matrix(Trajectories.Hem[,c("Spring_1", "Spring_2")]),
smoother='lowess',
trace=TRUE,
f = .7,
stretch=0)## Starting curve---distance^2: 45804778678
## Iteration 1---distance^2: 27732113
## Iteration 2---distance^2: 27728318
#The principal curve smoothed
Hem.pc.line <- as.data.frame(fit$s[order(fit$lambda),])
#Pseudotime score
Trajectories.Hem$PseudotimeScore <- fit$lambda/max(fit$lambda)if (cor(Trajectories.Hem$PseudotimeScore, Neurons.data@assays$SCT@data['Hmga2', Trajectories.Hem$Barcodes]) > 0) {
Trajectories.Hem$PseudotimeScore <- -(Trajectories.Hem$PseudotimeScore - max(Trajectories.Hem$PseudotimeScore))
}Trajectories.Pallial <- Neurons.data@meta.data %>%
select("Barcodes", "nUMI", "Spring_1", "Spring_2", "Lineage") %>%
filter(Lineage == "Pallial")fit <- principal_curve(as.matrix(Trajectories.Pallial[,c("Spring_1", "Spring_2")]),
smoother='lowess',
trace=TRUE,
f = .7,
stretch=0)## Starting curve---distance^2: 26984853690
## Iteration 1---distance^2: 22153700
## Iteration 2---distance^2: 22179462
## Iteration 3---distance^2: 22180297
#The principal curve smoothed
Pallial.pc.line <- as.data.frame(fit$s[order(fit$lambda),])
#Pseudotime score
Trajectories.Pallial$PseudotimeScore <- fit$lambda/max(fit$lambda)if (cor(Trajectories.Pallial$PseudotimeScore, Neurons.data@assays$SCT@data['Hmga2', Trajectories.Pallial$Barcodes]) > 0) {
Trajectories.Pallial$PseudotimeScore <- -(Trajectories.Pallial$PseudotimeScore - max(Trajectories.Pallial$PseudotimeScore))
}Trajectories.neurons <- rbind(Trajectories.Pallial, Trajectories.Hem)cols <- brewer.pal(n =11, name = "Spectral")
ggplot(Trajectories.neurons, aes(Spring_1, Spring_2)) +
geom_point(aes(color=PseudotimeScore), size=2, shape=16) +
scale_color_gradientn(colours=rev(cols), name='Speudotime score') +
geom_line(data=Pallial.pc.line, color="#026c9a", size=0.77) +
geom_line(data=Hem.pc.line, color="#cc391b", size=0.77)Neurons.data <- NormalizeData(Neurons.data, normalization.method = "LogNormalize", scale.factor = 10000, assay = "RNA")# Neurog2
p1 <- FeaturePlot(object = Neurons.data,
features = c("Neurog2"),
pt.size = 0.5,
cols = c("grey90", brewer.pal(9,"YlGnBu")),
reduction = "spring",
order = T) & NoAxes()
Trajectories.neurons$Neurog2 <- Neurons.data@assays$RNA@data["Neurog2", Trajectories.neurons$Barcodes]
p2 <- ggplot(Trajectories.neurons, aes(x= PseudotimeScore, y= Neurog2)) +
geom_point(aes(color= Lineage), size=0.5) +
scale_color_manual(values= c("#cc391b", "#026c9a")) +
geom_smooth(method="loess", n= 50, aes(color= Lineage)) +
ylim(0,NA)
# Tbr1
p3 <- FeaturePlot(object = Neurons.data ,
features = c("Tbr1"),
pt.size = 0.5,
cols = c("grey90", brewer.pal(9,"YlGnBu")),
reduction = "spring",
order = T) & NoAxes()
Trajectories.neurons$Tbr1 <- Neurons.data@assays$RNA@data["Tbr1", Trajectories.neurons$Barcodes]
p4 <- ggplot(Trajectories.neurons, aes(x= PseudotimeScore, y= Tbr1)) +
geom_point(aes(color= Lineage), size=0.5) +
scale_color_manual(values= c("#cc391b", "#026c9a")) +
geom_smooth(method="loess", n= 50, aes(color= Lineage)) +
ylim(0,NA)
# Mapt
p5 <- FeaturePlot(object = Neurons.data ,
features = c("Mapt"),
pt.size = 0.5,
cols = c("grey90", brewer.pal(9,"YlGnBu")),
reduction = "spring",
order = T) & NoAxes()
Trajectories.neurons$Mapt <- Neurons.data@assays$RNA@data["Mapt", Trajectories.neurons$Barcodes]
p6 <- ggplot(Trajectories.neurons, aes(x= PseudotimeScore, y= Mapt)) +
geom_point(aes(color= Lineage), size=0.5) +
scale_color_manual(values= c("#cc391b", "#026c9a")) +
geom_smooth(method="loess", n= 50, aes(color= Lineage)) +
ylim(0,NA)
p1 + p2 + p3 + p4 + p5 + p6 + patchwork::plot_layout(ncol = 2)Since we observe the first 25% of both trajectories are occupied by few, likely progenitor cells, we shift this cell along the axis
Pseudotime.intervals <- Trajectories.neurons%>%
select(Lineage, PseudotimeScore) %>%
mutate(Pseudotime.bins = cut(Trajectories.neurons$PseudotimeScore, seq(0, max(Trajectories.neurons$PseudotimeScore) + 0.05, 0.05), dig.lab = 2, right = FALSE)) %>%
group_by(Lineage, Pseudotime.bins) %>%
summarise(n=n())
ggplot(Pseudotime.intervals, aes(x=Pseudotime.bins, y=n, fill=Lineage)) +
geom_bar(stat = "identity", width = 0.90) +
theme(axis.text.x = element_text(angle = 45, hjust=1))+
scale_fill_manual(values= c("#cc391b", "#026c9a"))score <- sapply(Trajectories.neurons$PseudotimeScore,
FUN = function(x) if (x <= 0.2) {x= 0.2} else { x=x })
Trajectories.neurons$PseudotimeScore.shifted <- (score - min(score)) / (max(score) - min(score))# Neurog2
p1 <- FeaturePlot(object = Neurons.data ,
features = c("Neurog2"),
pt.size = 0.5,
cols = c("grey90", brewer.pal(9,"YlGnBu")),
reduction = "spring",
order = T) & NoAxes()
p2 <- ggplot(Trajectories.neurons, aes(x= PseudotimeScore.shifted, y= Neurog2)) +
geom_point(aes(color= Lineage), size=0.5) +
scale_color_manual(values= c("#cc391b", "#026c9a")) +
geom_smooth(method="loess", n= 50, aes(color= Lineage)) +
ylim(0,NA)
# Tbr1
p3 <- FeaturePlot(object = Neurons.data ,
features = c("Tbr1"),
pt.size = 0.5,
cols = c("grey90", brewer.pal(9,"YlGnBu")),
reduction = "spring",
order = T) & NoAxes()
p4 <- ggplot(Trajectories.neurons, aes(x= PseudotimeScore.shifted, y= Tbr1)) +
geom_point(aes(color= Lineage), size=0.5) +
scale_color_manual(values= c("#cc391b", "#026c9a")) +
geom_smooth(method="loess", n= 50, aes(color= Lineage)) +
ylim(0,NA)
# Mapt
p5 <- FeaturePlot(object = Neurons.data ,
features = c("Mapt"),
pt.size = 0.5,
cols = c("grey90", brewer.pal(9,"YlGnBu")),
reduction = "spring",
order = T) & NoAxes()
p6 <- ggplot(Trajectories.neurons, aes(x= PseudotimeScore.shifted, y= Mapt)) +
geom_point(aes(color= Lineage), size=0.5) +
scale_color_manual(values= c("#cc391b", "#026c9a")) +
geom_smooth(method="loess", n= 50, aes(color= Lineage)) +
ylim(0,NA)
p1 + p2 + p3 + p4 + p5 + p6 + patchwork::plot_layout(ncol = 2)ggplot(Trajectories.neurons, aes(x= PseudotimeScore.shifted, y= nUMI/10000)) +
geom_point(aes(color= Lineage), size=0.5) +
scale_color_manual(values= c("#cc391b", "#026c9a")) +
geom_smooth(method="loess", n= 50, aes(color= Lineage)) +
ylim(0,NA)rm(list = ls()[!ls() %in% c("Trajectories.neurons")])Progenitors.data <- readRDS("../ProgenitorsDiversity/Progenitors.RDS")table(Progenitors.data$Cell_ident)##
## Dorso-Medial_pallium Hem Medial_pallium
## 3451 1954 2719
To balance the number of progenitors in both domain we will only work with Hem and Medial_pallium annotated cells. Since we are using pallial cell to contrast CR specific trajectory we think this approximation will not significantly affect our analysis.
Progenitors.data <- subset(Progenitors.data, idents = c("Hem", "Medial_pallium"))p1 <- DimPlot(Progenitors.data,
reduction = "spring",
pt.size = 0.5,
cols = c("#e3c148", "#e46b6b")) + NoAxes()
p2 <- FeaturePlot(object = Progenitors.data,
features = "Angle.cc",
pt.size = 0.5,
cols = rev(colorRampPalette(brewer.pal(n =10, name = "Spectral"))(100)),
reduction = "spring",
order = T) & NoAxes()
p3 <- DimPlot(object = Progenitors.data,
group.by = "Revelio.phase",
pt.size = 0.5,
reduction = "spring",
cols = c(wes_palette("FantasticFox1")[1:3],"grey40",wes_palette("FantasticFox1")[5])) & NoAxes()
p1 + p2 + p3 + patchwork::plot_layout(ncol = 2)# Start with neurons data
Trajectories.all <- Trajectories.neurons %>% select(Barcodes, nUMI, Spring_1, Spring_2, Lineage)
Trajectories.all$Pseudotime <- Trajectories.neurons$PseudotimeScore.shifted + 1
Trajectories.all$Phase <- NA# Add progenitors data
Trajectories.progenitors <- Progenitors.data@meta.data %>%
select(Barcodes, nUMI, Spring_1, Spring_2) %>%
mutate(Lineage= ifelse(Progenitors.data$Cell_ident == "Medial_pallium", "Pallial", "Hem") ,
Pseudotime= Progenitors.data$Angle.cc,
Phase = Progenitors.data$Revelio.phase)Trajectories.all <- rbind(Trajectories.all, Trajectories.progenitors)
Trajectories.all$Phase <- factor(Trajectories.all$Phase, levels = c("G1.S", "S", "G2", "G2.M", "M.G1"))p1 <- ggplot(Trajectories.all, aes(Spring_1, Spring_2)) +
geom_point(aes(color=Pseudotime), size=0.5) +
scale_color_gradientn(colours=rev(brewer.pal(n =11, name = "Spectral")), name='Speudotime score')
p2 <- ggplot(Trajectories.all, aes(Spring_1, Spring_2)) +
geom_point(aes(color= Lineage), size=0.5) +
scale_color_manual(values= c("#cc391b", "#026c9a"))
p1 + p2p1 <- ggplot(Trajectories.all, aes(x= Pseudotime, y= nUMI/10000)) +
geom_point(aes(color= Phase), size=0.5) +
scale_color_manual(values= c(wes_palette("FantasticFox1")[1:3],"grey40",wes_palette("FantasticFox1")[5])) +
geom_smooth(method="loess", n= 50, fill="grey") +
ylim(0,NA)
p2 <- ggplot(Trajectories.all, aes(x= Pseudotime, y= nUMI/10000)) +
geom_point(aes(color= Lineage), size=0.5) +
scale_color_manual(values= c("#cc391b", "#026c9a")) +
geom_smooth(method="loess", n= 50, fill="grey") +
ylim(0,NA)
p1 + p2rm(list = ls()[!ls() %in% c("Trajectories.all")])Hem.data <- readRDS("../QC.filtered.cells.RDS")Neuro.trajectories <- CreateSeuratObject(counts = Hem.data@assays$RNA@data[, Trajectories.all$Barcodes],
meta.data = Trajectories.all)
spring <- as.matrix(Neuro.trajectories@meta.data %>% select("Spring_1", "Spring_2"))
Neuro.trajectories[["spring"]] <- CreateDimReducObject(embeddings = spring, key = "Spring_", assay = DefaultAssay(Neuro.trajectories))p1 <- FeaturePlot(object = Neuro.trajectories,
features = "Pseudotime",
pt.size = 1,
cols = rev(colorRampPalette(brewer.pal(n =11, name = "Spectral"))(100)),
reduction = "spring",
order = T) & NoAxes()
p2 <- DimPlot(object = Neuro.trajectories,
group.by = "Lineage",
pt.size = 1,
reduction = "spring",
cols = c("#cc391b", "#026c9a")) & NoAxes()
p3 <- DimPlot(object = Neuro.trajectories,
group.by = "Phase",
pt.size = 1,
reduction = "spring",
cols = c(wes_palette("FantasticFox1")[1:3],"grey40",wes_palette("FantasticFox1")[5])) & NoAxes()
p1 + p2 + p3rm(list = ls()[!ls() %in% c("Neuro.trajectories")])Neuro.trajectories<- NormalizeData(Neuro.trajectories, normalization.method = "LogNormalize", scale.factor = 10000, assay = "RNA")Neuro.trajectories <- FindVariableFeatures(Neuro.trajectories, selection.method = "disp", nfeatures = 3000, assay = "RNA")source("../Functions/functions_GeneTrends.R")
Plot.Genes.trend(Seurat.data= Neuro.trajectories,
group.by = "Lineage",
genes= c("Gas1","Sox2",
"Neurog2", "Btg2",
"Tbr1", "Mapt",
"Trp73", "Foxg1"))Plot.Genes.trend(Seurat.data= Neuro.trajectories,
group.by = "Lineage",
genes= c("Gmnc", "Mcidas",
"Foxj1", "Trp73",
"Lhx1", "Cdkn1a"))Plot.Genes.trend(Seurat.data= Neuro.trajectories,
group.by = "Lineage",
genes= c("Mki67", "Top2a",
"H2afx", "Cdkn1c"))# Transfer metadata
meta.data <- data.frame(Barcode= Neuro.trajectories$Barcodes,
Lineage= Neuro.trajectories$Lineage,
Pseudotime= Neuro.trajectories$Pseudotime,
Phase= Neuro.trajectories$Phase)
Annot.data <- new('AnnotatedDataFrame', data = meta.data)
# Transfer counts data
var.genes <- Neuro.trajectories[["RNA"]]@var.features
count.data = data.frame(gene_short_name = rownames(Neuro.trajectories[["RNA"]]@data[var.genes,]),
row.names = rownames(Neuro.trajectories[["RNA"]]@data[var.genes,]))
feature.data <- new('AnnotatedDataFrame', data = count.data)
# Create the CellDataSet object including variable genes only
gbm_cds <- newCellDataSet(Neuro.trajectories[["RNA"]]@counts[var.genes,],
phenoData = Annot.data,
featureData = feature.data,
lowerDetectionLimit = 0,
expressionFamily = negbinomial())gbm_cds <- estimateSizeFactors(gbm_cds)
gbm_cds <- estimateDispersions(gbm_cds)
gbm_cds <- detectGenes(gbm_cds, min_expr = 0.1)rm(list = ls()[!ls() %in% c("Neuro.trajectories", "gbm_cds", "Gene.Trend", "Plot.Genes.trend")])
gc()## used (Mb) gc trigger (Mb) max used (Mb)
## Ncells 3588903 191.7 6211198 331.8 6211198 331.8
## Vcells 102208537 779.8 404496316 3086.1 612662844 4674.3
pseudo.maturation.diff <- differentialGeneTest(gbm_cds[fData(gbm_cds)$num_cells_expressed > 80,],
fullModelFormulaStr = "~sm.ns(Pseudotime, df = 3)*Lineage",
reducedModelFormulaStr = "~sm.ns(Pseudotime, df = 3)",
cores = parallel::detectCores() - 2)# Filter genes based on FDR
pseudo.maturation.diff.filtered <- pseudo.maturation.diff %>% filter(qval < 1e-40)# Create a new pseudo-DV vector of 200 points
nPoints <- 200
new_data = list()
for (Lineage in unique(pData(gbm_cds)$Lineage)){
new_data[[length(new_data) + 1]] = data.frame(Pseudotime = seq(min(pData(gbm_cds)$Pseudotime), max(pData(gbm_cds)$Pseudotime), length.out = nPoints), Lineage=Lineage)
}
new_data = do.call(rbind, new_data)
# Smooth gene expression
Diff.curve_matrix <- genSmoothCurves(gbm_cds[as.character(pseudo.maturation.diff.filtered$gene_short_name),],
trend_formula = "~sm.ns(Pseudotime, df = 3)*Lineage",
relative_expr = TRUE,
new_data = new_data,
cores= parallel::detectCores() - 2)# Extract matrix containing smoothed curves for each lineages
Pal_curve_matrix <- Diff.curve_matrix[, 1:nPoints] #Pallial points
CR_curve_matrix <- Diff.curve_matrix[, (nPoints + 1):(2 * nPoints)] #CR points
# Direction of the comparison : postive ABCs <=> Upregulated in CR lineage
ABCs_res <- CR_curve_matrix - Pal_curve_matrix
# Average logFC between the 2 curves
ILR_res <- log2(CR_curve_matrix/ (Pal_curve_matrix + 0.1))
ABCs_res <- apply(ABCs_res, 1, function(x, nPoints) {
avg_delta_x <- (x[1:(nPoints - 1)] + x[2:(nPoints)])/2
step <- (100/(nPoints - 1))
res <- round(sum(avg_delta_x * step), 3)
return(res)},
nPoints = nPoints) # Compute the area below the curve
ABCs_res <- cbind(ABCs_res, ILR_res[,ncol(ILR_res)])
colnames(ABCs_res)<- c("ABCs", "Endpoint_ILR")
# Import ABC values into the DE test results table
pseudo.maturation.diff.filtered <- cbind(pseudo.maturation.diff.filtered[,1:4],
ABCs_res,
pseudo.maturation.diff.filtered[,5:6])# Extract Cajal-Retzius expressed genes
CR.res <- as.data.frame(pseudo.maturation.diff.filtered[pseudo.maturation.diff.filtered$ABCs > 0,])
CR.genes <- row.names(CR.res)
CR_curve_matrix <- CR_curve_matrix[CR.genes, ]## Cluster gene by expression profiles
Pseudotime.genes.clusters <- cluster::pam(as.dist((1 - cor(Matrix::t(CR_curve_matrix),method = "pearson"))), k= 5)
CR.Gene.dynamique <- data.frame(Gene= names(Pseudotime.genes.clusters$clustering),
Waves= Pseudotime.genes.clusters$clustering,
Gene.Clusters = Pseudotime.genes.clusters$clustering,
q.val = CR.res$qval,
ABCs= CR.res$ABCs
) %>% arrange(Gene.Clusters)
row.names(CR.Gene.dynamique) <- CR.Gene.dynamique$Gene
CR.Gene.dynamique$Gene.Clusters <- paste0("Clust.", CR.Gene.dynamique$Gene.Clusters)# Order the rows using seriation
dst <- as.dist((1-cor(scale(t(CR_curve_matrix)), method = "pearson")))
row.ser <- seriation::seriate(dst, method ="R2E") #"R2E" #TSP #"GW" "GW_ward"
gene.order <- rownames(CR_curve_matrix[get_order(row.ser),])
# Set annotation colors
pal <- wes_palette("Darjeeling1")
anno.colors <- list(lineage = c(Pallial_neurons="#026c9a", Cajal_Retzius="#cc391b"),
Gene.Clusters = c(Clust.1 =pal[1] , Clust.2=pal[2], Clust.3=pal[3], Clust.4=pal[4], Clust.5=pal[5]))
pheatmap::pheatmap(Diff.curve_matrix[gene.order,
c(200:1,#Pal
201:400)], #CR
scale = "row",
cluster_rows = F,
cluster_cols = F,
annotation_row = CR.Gene.dynamique %>% dplyr::select(Gene.Clusters),
annotation_col = data.frame(lineage = rep(c("Pallial_neurons","Cajal_Retzius"), each=200)),
annotation_colors = anno.colors,
show_colnames = F,
show_rownames = F,
fontsize_row = 8,
color = viridis::viridis(9),
breaks = seq(-2.5,2.5, length.out = 9),
main = "")anno.colors <- list(Cell.state = c(Cycling_RG="#046c9a", Differentiating_cells="#ebcb2e"),
Gene.Clusters = c(Clust.1 =pal[1] , Clust.2=pal[2], Clust.3=pal[3], Clust.4=pal[4], Clust.5=pal[5]))
col.anno <- data.frame(Cell.state = rep(c("Cycling_RG","Differentiating_cells"), each=100))
rownames(col.anno) <- 201:400
pheatmap::pheatmap(CR_curve_matrix[gene.order,],
scale = "row",
cluster_rows = F,
cluster_cols = F,
annotation_row = CR.Gene.dynamique %>% dplyr::select(Gene.Clusters),
annotation_col = col.anno,
annotation_colors = anno.colors,
gaps_col = 100,
show_colnames = F,
show_rownames = F,
fontsize_row = 8,
color = viridis::viridis(9),
breaks = seq(-2.5,2.5, length.out = 9),
main = "") ### Gene cluster trend
source("../Functions/functions_GeneClusterTrend.R")
Plot.clust.trends(Neuro.trajectories,
Lineage = "Hem",
Which.cluster = 1:5,
clust.list = Pseudotime.genes.clusters$clustering,
Smooth.method = "gam")CR.gostres <- gost(query = list("Clust.1" = CR.Gene.dynamique %>% filter(Gene.Clusters == "Clust.1") %>% pull(Gene) %>% as.character(),
"Clust.2" = CR.Gene.dynamique %>% filter(Gene.Clusters == "Clust.2") %>% pull(Gene) %>% as.character(),
"Clust.3" = CR.Gene.dynamique %>% filter(Gene.Clusters == "Clust.3") %>% pull(Gene) %>% as.character(),
"Clust.4" = CR.Gene.dynamique %>% filter(Gene.Clusters == "Clust.4") %>% pull(Gene) %>% as.character(),
"Clust.5" = CR.Gene.dynamique %>% filter(Gene.Clusters == "Clust.5") %>% pull(Gene) %>% as.character()),
organism = "mmusculus", ordered_query = F,
multi_query = F, significant = T, exclude_iea = T,
measure_underrepresentation = F, evcodes = T,
user_threshold = 0.05, correction_method = "fdr",
domain_scope = "annotated", custom_bg = NULL,
numeric_ns = "", sources = c("GO:MF", "GO:BP"), as_short_link = F)DNA_damage_GOterm <- CR.gostres$result[CR.gostres$result$term_id %in% c("GO:0008630", "GO:0030330", "GO:0031571", "GO:0006974", "GO:0006977","GO:0033554",
"GO:0044773", "GO:0042771", "GO:0042770", "GO:2001021", "GO:1902229"),]
DNA_damage_GOterm[,c(1,2,3,5,6,7,11)]## # A tibble: 8 × 7
## query significant p_value query_size intersection_size precision term_name
## <chr> <lgl> <dbl> <int> <int> <dbl> <chr>
## 1 Clust.1 TRUE 0.0172 93 3 0.0323 intrinsic …
## 2 Clust.2 TRUE 0.00622 182 28 0.154 cellular r…
## 3 Clust.4 TRUE 0.0159 121 4 0.0331 DNA damage…
## 4 Clust.4 TRUE 0.0283 121 4 0.0331 negative r…
## 5 Clust.4 TRUE 0.0288 121 21 0.174 cellular r…
## 6 Clust.4 TRUE 0.0339 121 3 0.0248 regulation…
## 7 Clust.4 TRUE 0.0367 121 5 0.0413 signal tra…
## 8 Clust.4 TRUE 0.0397 121 2 0.0165 DNA damage…
CR.gostres <- gost(query = as.character(CR.Gene.dynamique$Gene),
organism = "mmusculus", ordered_query = F,
multi_query = F, significant = T, exclude_iea = T,
measure_underrepresentation = F, evcodes = T,
user_threshold = 0.05, correction_method = "fdr",
domain_scope = "annotated", custom_bg = NULL,
numeric_ns = "", sources = c("GO:MF", "GO:BP"), as_short_link = F)DNA_damage_GOterm <- CR.gostres$result[CR.gostres$result$term_id %in% c("GO:0008630", "GO:0030330", "GO:0031571", "GO:0006974", "GO:0006977",
"GO:0044773", "GO:0042771", "GO:0042770", "GO:2001021", "GO:1902229"),]
DNA_damage_GOterm[,c(1,2,3,5,6,7,11)]## # A tibble: 10 × 7
## query significant p_value query_size intersection_size precision term_name
## <chr> <lgl> <dbl> <int> <int> <dbl> <chr>
## 1 query_1 TRUE 0.00107 510 11 0.0216 intrinsic…
## 2 query_1 TRUE 0.00430 510 7 0.0137 DNA damag…
## 3 query_1 TRUE 0.00505 510 4 0.00784 mitotic G…
## 4 query_1 TRUE 0.0118 510 30 0.0588 cellular …
## 5 query_1 TRUE 0.0118 510 3 0.00588 DNA damag…
## 6 query_1 TRUE 0.0184 510 6 0.0118 mitotic D…
## 7 query_1 TRUE 0.0207 510 5 0.00980 intrinsic…
## 8 query_1 TRUE 0.0365 510 9 0.0176 signal tr…
## 9 query_1 TRUE 0.0405 510 6 0.0118 negative …
## 10 query_1 TRUE 0.0496 510 4 0.00784 regulatio…
# Extract Pallial neurons trajectory genes
Pal.res <- as.data.frame(pseudo.maturation.diff.filtered[pseudo.maturation.diff.filtered$ABCs < 0,])
Pal.genes <- row.names(Pal.res)
Pal_curve_matrix <- Pal_curve_matrix[Pal.genes, ]## Cluster gene by expression profiles
Pseudotime.genes.clusters <- cluster::pam(as.dist((1 - cor(Matrix::t(Pal_curve_matrix),method = "spearman"))), k= 5)
Pal.Gene.dynamique <- data.frame(Gene= names(Pseudotime.genes.clusters$clustering),
Waves= Pseudotime.genes.clusters$clustering,
Gene.Clusters = Pseudotime.genes.clusters$clustering,
q.val = Pal.res$pval,
ABCs= Pal.res$ABCs
) %>% arrange(Gene.Clusters)
row.names(Pal.Gene.dynamique) <- Pal.Gene.dynamique$Gene
Pal.Gene.dynamique$Gene.Clusters <- paste0("Clust.", Pal.Gene.dynamique$Gene.Clusters)# Order the rows using seriation
dst <- as.dist((1-cor(scale(t(Pal_curve_matrix)), method = "spearman")))
row.ser <- seriation::seriate(dst, method ="R2E") #"R2E" #TSP #"GW" "GW_ward"
gene.order <- rownames(Pal_curve_matrix[get_order(row.ser),])
# Set annotation colors
pal <- wes_palette("Darjeeling1")
anno.colors <- list(lineage = c(Pallial_neurons="#026c9a", Cajal_Retzius="#cc391b"),
Gene.Clusters = c(Clust.1 =pal[1] , Clust.2=pal[2], Clust.3=pal[3], Clust.4=pal[4], Clust.5=pal[5]))
pheatmap::pheatmap(Diff.curve_matrix[gene.order,
c(200:1,#Pal
201:400)], #CR
scale = "row",
cluster_rows = F,
cluster_cols = F,
annotation_row = Pal.Gene.dynamique %>% dplyr::select(Gene.Clusters),
annotation_col = data.frame(lineage = rep(c("Pallial_neurons","Cajal_Retzius"), each=200)),
annotation_colors = anno.colors,
show_colnames = F,
show_rownames = F,
fontsize_row = 8,
color = viridis::viridis(9),
breaks = seq(-2.5,2.5, length.out = 9),
main = "")anno.colors <- list(Cell.state = c(Cycling_RG="#046c9a", Differentiating_cells="#ebcb2e"),
Gene.Clusters = c(Clust.1 =pal[1] , Clust.2=pal[2], Clust.3=pal[3], Clust.4=pal[4], Clust.5=pal[5]))
col.anno <- data.frame(Cell.state = rep(c("Differentiating_cells","Cycling_RG"), each=100))
rownames(col.anno) <- 200:1
pheatmap::pheatmap(Pal_curve_matrix[gene.order,200:1],
scale = "row",
cluster_rows = F,
cluster_cols = F,
annotation_row = Pal.Gene.dynamique %>% dplyr::select(Gene.Clusters),
annotation_col = col.anno,
annotation_colors = anno.colors,
gaps_col = 100,
show_colnames = F,
show_rownames = F,
fontsize_row = 8,
color = viridis::viridis(9),
breaks = seq(-2.5,2.5, length.out = 9),
main = "") ### Gene cluster trend
Plot.clust.trends(Neuro.trajectories,
Lineage = "Pallial",
Which.cluster = 1:5,
clust.list = Pseudotime.genes.clusters$clustering,
Smooth.method = "gam")#date
format(Sys.time(), "%d %B, %Y, %H,%M")## [1] "07 décembre, 2021, 18,44"
#Packages used
sessionInfo()## R version 4.1.2 (2021-11-01)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 20.04.3 LTS
##
## Matrix products: default
## BLAS: /usr/lib/x86_64-linux-gnu/blas/libblas.so.3.9.0
## LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.9.0
##
## locale:
## [1] LC_CTYPE=fr_FR.UTF-8 LC_NUMERIC=C
## [3] LC_TIME=fr_FR.UTF-8 LC_COLLATE=fr_FR.UTF-8
## [5] LC_MONETARY=fr_FR.UTF-8 LC_MESSAGES=fr_FR.UTF-8
## [7] LC_PAPER=fr_FR.UTF-8 LC_NAME=C
## [9] LC_ADDRESS=C LC_TELEPHONE=C
## [11] LC_MEASUREMENT=fr_FR.UTF-8 LC_IDENTIFICATION=C
##
## attached base packages:
## [1] splines stats4 stats graphics grDevices utils datasets
## [8] methods base
##
## other attached packages:
## [1] wesanderson_0.3.6 cowplot_1.1.1 ggExtra_0.9
## [4] RColorBrewer_1.1-2 dplyr_1.0.7 seriation_1.3.1
## [7] gprofiler2_0.2.1 monocle_2.22.0 DDRTree_0.1.5
## [10] irlba_2.3.3 VGAM_1.1-5 ggplot2_3.3.5
## [13] Biobase_2.54.0 BiocGenerics_0.40.0 Matrix_1.3-4
## [16] princurve_2.1.6 SeuratObject_4.0.4 Seurat_4.0.5
##
## loaded via a namespace (and not attached):
## [1] plyr_1.8.6 igraph_1.2.9 lazyeval_0.2.2
## [4] densityClust_0.3 listenv_0.8.0 scattermore_0.7
## [7] fastICA_1.2-3 digest_0.6.29 foreach_1.5.1
## [10] htmltools_0.5.2 viridis_0.6.2 fansi_0.5.0
## [13] magrittr_2.0.1 tensor_1.5 cluster_2.1.2
## [16] ROCR_1.0-11 limma_3.50.0 globals_0.14.0
## [19] matrixStats_0.61.0 docopt_0.7.1 spatstat.sparse_2.0-0
## [22] colorspace_2.0-2 ggrepel_0.9.1 xfun_0.28
## [25] RCurl_1.98-1.5 sparsesvd_0.2 crayon_1.4.2
## [28] jsonlite_1.7.2 spatstat.data_2.1-0 survival_3.2-13
## [31] zoo_1.8-9 iterators_1.0.13 glue_1.5.1
## [34] polyclip_1.10-0 registry_0.5-1 gtable_0.3.0
## [37] leiden_0.3.9 future.apply_1.8.1 abind_1.4-5
## [40] scales_1.1.1 pheatmap_1.0.12 DBI_1.1.1
## [43] miniUI_0.1.1.1 Rcpp_1.0.7 viridisLite_0.4.0
## [46] xtable_1.8-4 reticulate_1.22 spatstat.core_2.3-1
## [49] htmlwidgets_1.5.4 httr_1.4.2 FNN_1.1.3
## [52] ellipsis_0.3.2 ica_1.0-2 farver_2.1.0
## [55] pkgconfig_2.0.3 sass_0.4.0 uwot_0.1.10
## [58] deldir_1.0-6 utf8_1.2.2 labeling_0.4.2
## [61] tidyselect_1.1.1 rlang_0.4.12 reshape2_1.4.4
## [64] later_1.3.0 munsell_0.5.0 tools_4.1.2
## [67] cli_3.1.0 generics_0.1.1 ggridges_0.5.3
## [70] evaluate_0.14 stringr_1.4.0 fastmap_1.1.0
## [73] yaml_2.2.1 goftest_1.2-3 knitr_1.36
## [76] fitdistrplus_1.1-6 purrr_0.3.4 RANN_2.6.1
## [79] pbapply_1.5-0 future_1.23.0 nlme_3.1-153
## [82] mime_0.12 slam_0.1-49 rstudioapi_0.13
## [85] compiler_4.1.2 plotly_4.10.0 png_0.1-7
## [88] spatstat.utils_2.2-0 tibble_3.1.6 bslib_0.3.1
## [91] stringi_1.7.6 highr_0.9 lattice_0.20-45
## [94] HSMMSingleCell_1.14.0 vctrs_0.3.8 pillar_1.6.4
## [97] lifecycle_1.0.1 spatstat.geom_2.3-0 combinat_0.0-8
## [100] lmtest_0.9-39 jquerylib_0.1.4 RcppAnnoy_0.0.19
## [103] bitops_1.0-7 data.table_1.14.2 httpuv_1.6.3
## [106] patchwork_1.1.1 R6_2.5.1 promises_1.2.0.1
## [109] TSP_1.1-11 KernSmooth_2.23-20 gridExtra_2.3
## [112] parallelly_1.29.0 codetools_0.2-18 MASS_7.3-54
## [115] assertthat_0.2.1 withr_2.4.3 qlcMatrix_0.9.7
## [118] sctransform_0.3.2 mgcv_1.8-38 parallel_4.1.2
## [121] grid_4.1.2 rpart_4.1-15 tidyr_1.1.4
## [124] rmarkdown_2.11 Rtsne_0.15 shiny_1.7.1
Institute of Psychiatry and Neuroscience of Paris, INSERM U1266, 75014, Paris, France, matthieu.moreau@inserm.fr↩︎